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We study homoclinic orbits of the Swift-Hohenberg equation near a 
Hamiltonian-Hopf bifurcation. It is well known that in this case the nor- 
mal form of the equation is integrable at all orders. Therefore the dif- 
ference between the stable and unstable manifolds is exponentially small 
and the study requires a method capable to detect phenomena beyond 
all algebraic orders provided by the normal form theory. We establish 
an asymptotic expansion for a homoclinic invariant which quantitatively 
describes the transversality of the invariant manifolds. We perform high- 
precision numerical experiments to support validity of the asymptotic ex- 
pansion and evaluate a Stokes constant numerically using two independent 
methods. 



1 The generalized Swift-Hohenberg equation 

The generalized Swift-Hohenberg equation (GSHE), 



is widely used to model nonlinear phenomena in various areas of modern Physics 
including hydrodynamics, pattern formation and nonlinear optics (e.g. [11 110|V 
This equation (with k = 0) was originally introduced by Swift and Hohen- 
berg [18] in a study of thermal fluctuations in a convective instability. 

In the following we consider u to be one dimensional and study stationary 
solutions of ([!]) which satisfy the ordinary differential equation 
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Figure 1: Two primary symmetric iiomoclinic solutions of the scalar stationary GSHE 
(e= -0.05). 

Obviously this equation has a reversible symmetry (if u{x) satisfy the equation 
then u{—x) also does). It is well known that for small negative e this equation 
has two symmetric homoclinic solutions [9 similar to the ones shown on Fig- 
ure [1] In this paper we study transversal! ty of the homoclinic solutions, which 
implies existence of multi-pulse homoclinic solutions and a small scale chaos. 

In order to describe the homoclinic phenomena it is convenient to rewrite 
the equation ^ in the form of an equivalent Hamiltonian system [3 [T3] : 



91 = 92 



Pi =P2- £91 - Kql + ql 



(3) 



P2 = -Pl , 



where the variables are defined by the following equalities 



u 91, 



u = 92, 



{u' + u'") = Pi and u + u" = p2 



(4) 



and the Hamiltonian function has the form 




(5) 
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The system ^ is reversible with respect to the involution, 



S ■■ {qi,q2,Pi,P2) {qi,-q2,-pi,P2)- 

The origin is an equiblibrium of the system and the eigenvalues of the linearized 
vector field are 



If e < 0, the eigenvalues form a quadruple ±/3e ± ia^ where 




At e = the eigenvalues collide forming two purely imaginary eigenvalues ±i 
of multiplicity two. Moreover, the corresponding linearization of the vector 
field is not semisimple. Thus, the equilibrium point of system ([3]) undergoes 
a Hamiltonian-Hopf bifurcation described in the book [TB] (see also [TT]). In 
general position there are two possible scenarios of the bifurcation depending 
on the sign of a certain coefficient of a normal form. In the Swift-Hohenberg 
equation both scenarios are possible and depend on the value of the parameter 
K. In this paper we will consider the case when the equilibrium is stable at the 
moment of the bifurcation (see [HI [15] for more details) which corresponds to 

|k| > -^/ll as shown in [T]. 

When e < is small, the equilibrium is a saddle-focus and the Stable Mani- 
fold Theorem implies the existence of two-dimensional stable Wf and unstable 
W" manifolds for the equilibrium point. These manifolds are contained inside 
the zero energy level of the Hamiltonian i/j . The normal form is integrable and 
the corresponding stable and unstable manifolds coincide (see also discussion 
in [12] for the reversible set up). It can be also deduced that the Hamiltonian 
system has two symmetric homoclinic orbits [I] but no conclusion about their 
transversality can be made using the normal form theory. In this paper we will 
study the splitting of the stable and unstable manifolds which happens beyond 
all orders of the normal form theory. We find a natural way to normalise vectors 
tangent to Wf and W" at the primary symmetric homoclinic point. The 
main result of this paper is the following asymptotic formula for the value of 
the standard symplectic form on this pair of vectors: 

n{v'',v') = e" W (cjo(k) + 0(e)) . (6) 

This formula implies the transversality of the homoclinic orbit for all small values 
of e provided the splitting coefficient wo(k) does not vanish. This constant is 
related to a purely imaginary Stokes constant, and Figure [D gives an idea about 
its behaviour as a function of the parameter k. 

We note that the derivation of the equation (|6]) does not rely on the specific 
form of the Swift-Hohenberg equation and the same asymptotic expression is 



3 




valid near the Hamiltonian-Hopf bifurcation in a generic analytic family. The 
value of the Stokes constant comes from the study of the Hamiltonian at the 
exact moment of the bifurcation. 

The original Hamiltonian ([5]) can be seen as a perturbation of an integrable 
Hamiltonian H" which can be derived from the order n normal form theory. 
Glebsky and Lerman 9 proved the persistence of the reversible homoclinic 
orbits for the original system ([3]) when e < is small. Their proof is based 
on an application of the implicit function theorem. As a matter of fact, this 
result follows from a more general study concerning a 1:1 resonance in four 
dimensional reversible vector fields (see [E])- Also the paper |^ conjectures 
that the stable and unstable manifolds should intersect transversely yielding, in 
particular, the existence of countably many reversible homoclinic orbits. These 
orbits are known as multisolitons for the Swift-Hohenberg equation and they 
have been the subject of study in several works (see [5] and the references 
therein) . 

In this paper we study the non-degenerate scenario of the bifurcation and 
note that the degenerate case leads to some interesting phenomena including 
"homoclinic snaking" [501 US] ■ 

The rest of the paper is organised in the following way. In Section 11.11 we 
discuss the definition of a homoclinic invariant which provides a very convenient 
tool for measuring the splitting of invariant manifolds. In section [1.21 we review 
some facts from the normal form theory which will be necessary for the expo- 
sition of our results. Section 11.41 introduces the main result of this paper: an 
asymptotic formula which describes the splitting of invariant manifolds of the 
stationary Swift-Hohenberg equation near the Hamiltonian-Hopf bifurcation. In 
particular this formula implies existence of two transversal primary homoclinic 
orbits and gives a sharp lower (exponentially small) bound for the splitting of 
the invariant manifolds. The derivation of the asymptotic formula is not com- 
pletely rigorous and we perform a set of high-precision numerical experiments 
in order to confirm its validity. Moreover, similar to many other problems which 
involve exponentially small splitting of invariant manifolds [7 , the asymptotic 
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formula contains a "splitting" constant which conies from an auxiliary problem 
and requires numerical evaluation. The results of our numerical experiments 
are reported in Sections [2] and [3] 

1.1 Homoclinic invariant 

In a study of homoclinic trajectories, both numerical and analytical, it is usu- 
ally important to have a convenient basis in the tangent space to the stable 
and unstable manifolds. Below we provide a definition adapted to our problem. 
This definition can be of independent interest as it can be easily extended onto 
hyperbolic equilibria of higher dimensional systems (not necessarily Hamilto- 
nian) . 

Suppose that the origin is an equilibrium of a Hamiltonian vector field Xh 
and that ±/3 ± ia are the eigenvalues of DXh{Q)- Then the origin has a two 
dimensional stable manifold. According to Hartman 11; the restriction of the 
vector field on Wi^^ can be linearised by a change of variables. In the polar 
coordinates the linearised dynamics on Wf^^ take the form: 

f = —/3r if — a . 

It is convenient to introduce z = — In r so that 

z = /3. 

Then the local stable manifold is the image of a function 

r : {{ifi, z):ipeS\z>- logro} ^ M"^ 

where rp is the radius of the linearisation domain and is the unit circle. Since 
maps trajectories into trajectories we can propagate it uniquely along the 
trajectories of the Hamiltonian system using the property 

r'{ip + at,z + f3t) = <i>*Hor''{ip,z) (7) 

where <I>|f is the fiow defined by the Hamiltonian equation. Note that 

r'*(^ + 27r,z) = r%ip,z) 

since ip is the angle component of the polar coordinates. Moreover, 

lim r'*(^,z) =0. 

Differentiating along a trajectory we see that it satisfies the non-linear PDE: 

ad^r + (3d,T = XniT) . (8) 

Each of the derivatives dzT'^ and dipT" defines a vector field on W. The equation 
([7} implies that dzT'' and dipV are invariant under the restriction of the fiow 
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We can define F" applying the same arguments to the Hamiltonian —H. In 
this case it is convenient to set z = In r to ensure that F" satisfies the same PDE 
as F*. In a reversible system with a reversing involution 5*, it is convenient to 
set 

F"((^,z) ^ SoV'{-ip,-z). 

Now suppose that the system has a homoclinic trajectory 7/1. Let us choose a 
point pft, e 7/,.. The freedom in the definition allows us to assume that p/j = 
F'*(0, 0) = F"(0, 0) without loosing in generality. This condition completely 
eliminates the freedom from the definition of F" and . 

In a Hamiltonian system the symplectic form provides a natural tool for 
studying transversality of invariant manifolds. The homoclinic invariant uj is 
defined by 

t^ = l}(a^F"(0,0),5^F^(0,0)). (9) 

It is relatively straightforward to check that w is an invariant: the defini- 
tion leads to the same value for all points of the homoclinic trajectory = 
{$^(pft) : t S M}. We also note that the definition of w does not depend 
on the choice of coordinates. Moreover, since F** and belong to the energy 
level H — i?(0), which is three-dimensional, the inequality a; 7^ implies the 
transversality of the homoclinic trajectory. 

We note that we defined two vectors tangent to and another two vectors 
tangent to at \ih € n . So we could use 

:=r!(aj"(0,0),a,F^(0,0)) 

instead of lj. But this invariant is not independent. Indeed, using the equality 

a5^F"(0, 0) + PdzV'iQ, 0) = ad^T'iQ, 0) -I- /35,F'(0, 0) 

(both are equal to XHiVh)) and the fact that W"^'^ are Lagrangian manifolds, 
we can easily deduce 

a^uj + = . 

Wc note that our definition of the homoclinic invariant is a natural extension 
of the invariant defined for homoclinic orbits of area-preserving maps [7] . 

In the case of the Swift-Hohenberg equation the system of PDE ([8]) can be 
conveniently replaced by a single scalar PDE of higher order obtained from ^ 
by replacing dx with the differential operator 

d = a^d^ -I- P^d^. 

Let us use uf- to denote the first component of and respectively, then 
satisfies the equation 

{l + d'^fu = eu + Ku'^-u^. (10) 

Its other components can be restored using The Swift-Hohenberg equation 
is reversible and we assume that 
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and rg(0, 0) = r"(0, 0) G S is the primary symmetric homoclinic point. Then 
the formula for the homochnic invariant can be rewritten in terms of u~: 



where the derivatives are evaluated at (tp, z) = (0, 0). 



(11) 



1.2 Normal form of the Swift-Hohenberg equation 

The most convenient description of the bifurcation is obtained with the help of 
the normal form. As a first step the quadratic part of the Hamiltonian ([S]) is 
normalised with the help of a linear symplectic transformation (similar to |3]): 



T = 



I 

1/4 %/2 

V 



-1/4V2 



-V2 



-1/2 






\ 

1/2 V2 


) 



which transforms ([5]) into 



Hc = - - qiP2) 



1 

16 



V2 
16 



1. 2 
2(91- 

96 



92; 



V2 3 



16 



Pt 



^'72Pie- 



32 



2 2 
92P1 



a/2 2 

— K(72Pi + 
1 3 1 4 

32'^^^^ - 256'^^ 
(12) 



where we keep the same notation for the variables. Note that the involution S 
in the new coordinates takes the form 



S ■■ iqi,q2,Pi,P2) {-qi,q2,pi,-p2)- 



(13) 



Now, with the quadratic part in normal form, we can apply the standard normal 
form procedure to normalize the Hamiltonian (|12p up to any order: There is 
a near identity canonical change of variables which normalizes all terms of 
order less than equal to n and transforms the Hamiltonian to the following form: 



higher order terms 



(14) 



where 



with 



h = q2Pi - qiP2, 



-I2 + 



3i+2j"+2;>4 
t+j>l 



ql 



ql 



P\ 



pI 
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This normalization preserves the reversibihty with respect to the involution (IT^ . 
In the case of the GSHE the normal form up to the order five has the form (see 
Appendix 1^ for more details about the change of variables) 

The leading part of the normal form includes two parameters which can be 
explicitly expressed in terms of the original parameter k: 

. 19 2 3 \ , „ / 65 2 3 

71 = 4 K.^ and u = 2[ 

' V576 1287 V864 64 

The geometry of the invariant manifolds depends on the sign of rj 16J. In the 
case of GSHE, if 

[27 

then ?7 > [9] , and the truncated normal form has a continuum of homoclinic 
orbits among which exactly two are reversible, i.e., symmetric with respect to 
the involution ([13]). 

In order to describe the geometry of the invariant manifolds near the bifur- 
cation it is convenient to introduce the new parameter e = — 4i5^ and perform 
the standard scaling: 

qi^S'^Qi, q2 = S'^Q2, pi^SPi, P2 = 5P2. 

This change of variables is not symplectic, nevertheless it preserves the form of 
the Hamitltonian equations since the symplectic form gains a constant factor 
5'^, so we have to multiply the Hamiltonian by in order to return back to 
the standard symplectic form. The Hamiltonian i/" is transformed into. 



This Hamiltonian system has an equilibrium at the origin characterized by a 
quadruple of complex eigenvalues ±za„_e ± (3n,e, where a„.e = 1 + ^5^ + 0{5'^) 

The equilibrium has a two dimensional stable and two dimensional unstable 
manifolds. It is convenient to parametrize these manifolds by solutions of the 
partial differential equation: 

(a„,,a^ + Y„ = X^j(T„). (15) 

The function T„(iy9, z) is real-analytic, converges to zero as z — )■ ±00 and is 
27r-periodic in (p. Taking into account the rotational symmetry of the normal 
form Hamiltonian, we can look for the solution of this equation in the form: 

Tn{(p,z) = (i?„(z)cos(6'„((^, z)),i?„(z)sin(6'„((^, z)), 
r„(z) cos(6'„((/3, z)), r„(z) sin(6'„((^, z))) 
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where Rn{z), rn{z) and 9n{'f,z) are real analytic functions. In particular, for 
n = 5 we get the following system of equations: 

/35.i?^ = -<5r5 (1 - wl) , /35,er^ = -5R^ , 

(a5,.a^ + /35.e9.) 05 = 1 + y (1 - t^rl) . 



From these equations we conclude, if 



5 



2 



then 



, _ 1 sinh z 

r5 = \ 1 — , ^5 



T] cosh z ' V 7? cosh^ z 

^^t^ 2 , sinhz 



2 J T] coshz 

We see that (r(z),i?(z)) runs over a homoclinic loop when z varies from — oo 

to +00. 

In general the parameterization T„ is the unique solution of (1151) such that 
-Rn(O) = and 6n{ip,0) = ip. Thus, T„((/3, z) belongs to the symmetry plane 
associated with the involution ()13p if and only if z = and (p = or ip = it. 
Therefore, there are exactly 2 symmetric homoclinic points. Let us call these 
homoclinic orbits the primary reversible homoclinic orbit. 

1.3 Stokes constant 

Let us study invariant manifolds of GSHE for e = 0. It is convenient to 
parametrise these invariant manifolds by solutions of the following PDE: 

(1 + {d^ + drf)^U = KU^ - . (16) 

It can be shown that this equation has an analytic solution u = Uq with the 
following asymptotic behaviour: 

u^i^,r)^^ + ^ + Oir-^) (17) 

in the set 

rep-g^ ={t : |arg(T + r)| >0o} , 
where is a small fixed constant and r is sufficiently large and 

The function Uq is 27r-periodic in ip. 
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The equation (jl6p has a second sohition u = Uq with 

It has the same asymptotic behaviour as Uq but is defined in a different sector, 
more precisely, it is defined for r such that — r G T^rBo' "^^^ solutions have 
a common asymptotics on the intersection of their domains but they do not 
typically coincide. The difference of these two solutions can be described in the 
following way. We can restore 4-dimensional vectors using equations (|4]) 
with ' replaced by + dr- In particular, the first component of coincides 
with . The functions are parameterizations of the stable and unstable 
manifolds and satisfy the following non-linear partial differential equation, 

{d^+dr)T = XHA^) (19) 

Let 

Ao(¥',t) = r+(^,T)-ro((^,T) 

and 

0o(¥',r)=l^(Ao(^,T),9^r+((^,r)), 

where $7 is the standard symplectic form. 

It can be shown that there is a constant 0o(k) such that 

eo{ip, t) = eo(K)e-'("~'^) + O(e-(2-^0)^(r-p)-) (20) 

as Imr — > — oo and for very small eq > 0. The constant Oo{k,) is known as 
the Stokes (or splitting) constant. The Stokes constant can be defined by the 
following limit: 

eo(^):= lim 0o{ip,T)e''-^-'^K (21) 

Im(T)— ^ — oo 

We note that the value of the Stokes constant cannot be obtained from our 
arguments. Fortunately the numerical evaluation of this constant is reasonably 

easy. Figure [2] shows the values of Im0o(K) plotted against k for k > kq = 

The picture suggests that the Stokes constant vanishes infinitely many times and 

that its zeros accumulate to kq- 



1.4 Asymptotic formula for the homoclinic invariant 

In this section we derive an asymptotic formula for the homoclinic invariant of 
the primary symmetric homoclinic orbit using the complex matching method. 

In the previous section we obtained an approximation of the separatrix in 
the normal form coordinates. Transforming T^{ip,z) back to the original coor- 
dinates we obtain the following approximation: 

, s. 1 COS {if) ^ . 

^ / 9k-Hkcos(2(p) 1 sin((^) sinh(z) \ ^ ^^^3^ 



I877 cosh (z) ^/rj \tl 2/ cosh (z) 
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where e — — 4^^. Since the function in the right- hand-side of the equation 
is even, it also approximates the stable separatrix represented by uf{ip,z) = 
u~ {—if, —z). A more accurate approximation with a 0((5") error can be obtained 
with the help of higher order normal form theory, but naturally none of those 
approximations can distinguish between the stable and unstable separatrices 
and we come to the conclusion that 

u-iip,z)~ut{ip,z)^Oin 

for all n. Of course the constant in this upper bound may depend on the point 
{(p,z). Therefore, the difference between the stable and unstable parametrisa- 
tion cannot be detected using power series of the perturbation theory, and we 
say it is beyond all algebraic orders. A rather standard approach to the prob- 
lem is based on studying the analytical continuation of the parametrisations and 
looking for places in the complexified variables where the leading orders of the 
approximation (j22l) grow significantly. We note that the variables z and ip play 
different roles, in particular we assume that ip is kept real or, more precisely, in 
a fixed narrow strip around the real axis. 

It is easy to see that the leading orders of u~ have poles at z = + kin for 
any integer k. In the following we study the behaviour of the parametrisations 
near the singular point z — i^. The first step is to re-expand the functions in 
Laurent series around the singularity and introduce a new variable 

Substituting this new variable into (|22p and expanding around r = we con- 
clude that 

u:i^, + ^f) . (Me) + ^ + o(,-3)^ + oie) (24) 

where Pi and P2 are the same as in and the error terms come from the 
analysis of the next order corrections. In this analysis we consider the terms in 
([22]) which are most divergent and in this way obtain the essential behaviour of 
u~ around the singularity. 

Transforming the equation ^TU\\ to the variable , setting e = and noting 
that ao = 1), we obtain equation (jl6p considered in the previous subsection. 
The following method is known as "complex matching" and is based on the 
observation that u J approximate uf- in a region where \z — ij\ is small but r 
is still large. Taking into account (|24p we conclude that 

^Tif^^J + = u^if^T) +0{e) , 

in a neighbourhood of a segment of the imaginary axis where Sr is large nega- 
tive. In a rigorous justification of the method we use the interval — i?loge~^ < 
Qr < —R, where _R is a large constant. 
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In order to derive an asymptotic formula for the homoclinic invariant, we 
consider an auxiliary function defined by 

where A((^, z) = r"(iy9, z) — {ip, z) and Vl is the standard symplectic form. The 
homoclinic invariant of the primary homoclinic orbit is defined by ^ which 
takes the form 

w = i](9^r^(o,o),a^r:(o,o)). (25) 

Differentiating the definition of 6 at the origin and taking into account that 
A(0, 0) = we get the relation: 

w = a^e(o,o). 

We claim that in a neighbourhood of the point t ~ — ilog(e^^) the following 
estimate holds: 

A(^,^r + zf)--Ao(^,T) + 0(62) (26) 

which leads to 

e(^, z) = -0o{p, t) + 0(e2) = -e-^(-'P)eo{K) + 0{e^) . (27) 

Let d — a^d^ + Pedz- It is easy to see that if vi and V2 both satisfy the first 
order variational equation 

dv = DXH,(Tl)v, (28) 

then 

9fl(vi, V2) = fi(9vi, V2) + ri(vi, 9V2) 

= n{DXHSTi)wi,w2) + n{wi,DXHSK)^2) = o. 

Therefore r2(vi((y9, z), V2((y5, z)) can be considered as a function of a single vari- 
able f{acZ — Pe(p). Moreover if vi^2 are 27r-periodic in ip, the function / is 
27r^-periodic in z. 

We note that d^pT^ satisfies the variational equation ((28)) . and A is an ap- 
proximate solution with the error of the order of 0(|A((^,z)p). Taking into 
account that the splitting of separatrices is rather small, we continue our argu- 
ments neglecting this error. Then there is a function / such that 

Qi'PiZ) = f{a^z - P^ip) 

inside the domain of @, which implies that / can be extended by periodicity 
onto the strip |5(z)| < f — We expand the function / into Fourier series, 
i.e., 

e(^,z) = 5]/,e"^(t^--). 

fcez 
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The coefficients of tfie series can be expressed in terms of Fourier integrals: 

271- ,9, 



/fc = ^ / e *V^e(0,z)dz. (29) 



Following the common procedure of Fourier Analysis, we shift the contour of 
integration to 3z = J — ^loge^^, change the variable to (|23| and use the 
estimate ([27]) to get 

/-I = -e"2/Jc (eo(«;)+0(e)) , (30) 

fi — f-i and there is a positive constant C such that 

|/fe| < Ce'^l'^le"' '2/3, for|fc|>2. 

Substituting these estimates into the Fourier series we get that for real values 
of ifi, z 

e(v3, z) = -2e" w |eo| cos ( - - arg(eo) ) + 0(e" We) , 

): { ... 

d^e{ip, z) = -2e- W lOol sin f ^z - - arg(eo) 1 + Oie'^e) . 

Since 8(0,0) = for all e then arg(9o) = ±f , i.e., the Stokes constant is a 
purely imaginary number and equation ^ follows directly. 

We note that the integrability of the normal form allows us to repeat the 
arguments with more accurate approximations of the separatrices, the result of 
this consideration leads to the conjecture that 

u;{e)-e~W7Y^u;ke' (32) 



fc>0 



where ujq = 2Im(8o(K))- 



2 Computation of the Stokes constant 

Since the arguments involved in the derivation of the asymptotic formulas are 
not completely rigorous, we have developed numerical methods to check the 
validity of our results. The procedure is based on comparison of two different 
methods for evaluation of the Stokes constants. The first method relies on the 
definition (|2ip and involves the GSHE with e = only. The second method 
evaluates the homoclinic invariant for e ^ and relies on the validity of the 
asymptotic expansion (|32|) to extrapolate the values of the (normalised) homo- 
clinic invariant towards e = in order to get luq. 
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Figure 3: Graph of log ( 73^;'^^)^' ) 

2.1 A method for the computation of the Stokes constant 

Let us describe the first method for computing the Stokes constant. We set 
T = —ia, if = and rewrite equation (1201) in the form: 

eo = 0o(O,-^^T)e" + o(e-(l-^'')") (33) 

Then we proceed as follows. 

(i) The first step is to construct a good approximation of stable and unstable 
manifolds. This approximation is given by a finite sum of the unique 
formal separatrix Tq of p9)) . 

N 

fe=l 

where 

fc 

Tkiv) = rfeje^'^-^ with Tfe,, e 

Note that Fjv approximates the parameterizations in the following 
sense 

The natural number N can be chosen using the astronomers recipe. It 
simply chooses N such that for fixed r and ip it minimizes |rfc((/?)T~^~^ | , 
that is, the least term of the series ro(i^,T) (see Figure [3]). 

(ii) A point on the unstable manifold (resp. stable manifold) can be repre- 
sented in the coordinates t). In order to obtain a point close to the 
unstable manifold we fix a positive real number a G IR+ and a sufficiently 
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large d € K"*" and define Zq = TN{—d,—ia — d) and a tangent vector 
Vq = d^TN{^d, —i<7 — d). Analogously, for the stable manifold we define 
z'^ ~ TN{d, —ia + d) and Vq = dipr]^{d, —ia + d). 

(iii) The next step is to measure the difi'erence of stable and unstable manifold 
at the point (</?, t) = (0,— iu). Taking into account the periodicity in Lp 
we set d equal to a multiple to 27r and integrate numerically the system, 

z' ^ Xho (z) (oa^ 

v'^DXh,{z)v ^"^^^ 



forward in time with t € [0, d] and initial conditions z (0) = Zq , w (0) = 
Vq and then backward in time with t d [— d, 0] and initial conditions 

z+iO) = z+,v+{0) = v+. 

(iv) Finally we evaluate, 

e{a) = n(z+{-d) - z"(d),w-(d))e'" (35) 



Remark 1 The stable and unstable manifolds have the same asymptotic expan- 
sion, hence the difference z~^{—d) — z~{d) is exponentially small, i.e. comparable 
with e" . Thus the system p4p has to be integrated with great accuracy. In the 
case of GSHE an excellent integrator can be constructed using a high order Tay- 
lor series method. 



2.2 Numerical results 

In all current computations we have used a Taylor series method, which is 
incorporated in the Maple Software, to integrate the equations of motion ((34|) . 
The method uses an adaptive step procedure controlled by a local error tolerance 
which was set to lO"'^, where D is the number of significant digits used in the 
computations. The order of the method has been automatically defined using 
the formula max(22, [1.5Z5J). Having fixed k = 2 we have computed the first 
45 coefficients of the formal separatrix Tq with 60 digits precision. The error 
committed by the approximation r^r is approximately of the order of the first 
missing term. Using double precision (16 digits) we have integrated numerically 
the equations ([M]) to obtain Q{a) for values of a uniformly distributed in the 
interval [20,28.89]. The initial conditions were computed using d = 3507r and 
the first 9 terms of the formal series Tq. The results are depicted in Figure 21 
The expected errors are bounded by the red curves. This implies in particular 
that the method is numerically stable, that is, the propagation errors due to 
integration do not increase drastically. There are several sources of errors that 
affect the accuracy of the computation of the Stokes constant, namely: 

• Approximation of stable and unstable manifolds given by the function Fjv; 

• Errors due to the numerical integration; 
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Figure 4: The top figure represents the graph of the function lm(0(a))e'^ — 
10.472161956944 and the bottom figure represents the graph of the function 
Re(0((j))e'^. When a is around 25 the rounding errors become visible and the conver- 
gence stops. The red curves represent the magnitude of the rounding errors. 
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D 




Rc(e(<Tn) 


lm(0{^-^ 


)) 


16 


24.68 


2.7e-05 


10.4721f 


5143901571 


20 


29.46 


7.8e-07 


10.472161953423286113 


24 


34.21 


1.6C.-08 


10.472ie 


519569069446924024 


28 


38.95 


3.1e-10 


10.4721f 


5195694413924682820786 


32 


43.67 


5.3e.-12 


10.4721f 


J1956944396725504278408504 


36 


48.37 


8.5e-14 


10.4721f 


519569443983419527788851129556 


40 


53.07 


1.2e.-15 


10.4721t 


J195694439S35S 129892633 11456886391 


44 


57.76 


1.8e-17 


10.4721f 


5195694439S358284 1806844684678196221 91 


48 


62.45 


2.6e.-19 


10.4721f 


J1956944398358285508435672590071 7201861670 


52 


67.12 


3.5e-21 


10.4721f 


5195694439835828552130242825730920048239485015 


56 


71.80 


4.7e.-23 


10.4721t 


51956944398358285521430879142372532568396894067732 


60 


76.46 


6.2C-25 


10.4721f 


519569443983582855214320209319731283197852962601326570 


64 


81.13 


8.0C-27 


10.4721f 


5195694439835828552143203166495538939445255794702026972749 


68 


85.79 


l.Oc-28 


10.4721f 


J195694439835828552 143203 1900047829633854060398152634432422925 



Table 1: Stokes constant evaluated at the optimum a* for dilFcrcnt computer 
precisions. In the computations we have used d ~ SSOvr and = 40 



Rounding errors. 



The first and the second source of errors can be made small compared to the 
rounding errors, which can be roughly estimated by, 



C 



■10 



(36) 



where D is the number of digits used in the computations and C is a real positive 
constant which reflects the propagation of rounding errors. Using this estimate 
we have provided bounds for the rounding errors which can be observed in 
Figure 131 The constant C can be estimated by fitting the function ([55)1 to the 



points 



e(a) 



for CT > 25. Using the method of least squares we have concluded 
that C is approximately 16.7. 

With double arithmetic precision the method previously described allows 
the computation of 7 to 8 correct digits of the Stokes constant Qq. In fact the 
rounding errors in computing G)((t) from formula (j35p grow accordingly to (j36l) 
whereas the neglected terms of the formula ((33|) decrease like Cie^"' , where Ci is 
some positive constant. Hence the optimal is attained when both contributions 
are of the same order. The constant Ci can be estimated by fitting the function 



Co + Cie to the points 



e(a) 



for a < 24. Using the method of least squares 
we have obtained that Ci is approximately 17305.75. Using this information we 
can determine the value a* where both contributions are essentially of the same 
order. This means that a* must satisfy the equation. 



a\2 



C 



\v\a'C, 



-10-^ 



which implies. 



eo-e(a*) 



816 



-IQ- 



In this way it is possible to obtain 8 correct digits for the Stokes constant 
using only double precision. In Table [T] we have listed the values of 9(cr*) 
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d\N 


10 


20 


30 


IOOtt 


10.47216215179386 


10.47216215183208 


10.47216215181955 


ISOvr 


10.47216131335742 


10.47216131335746 


10.47216131335772 


2007r 


10.47216144775669 


10.47216144775671 


10.47216144775682 


2507r 


10.47216149546998 


10.47216149546998 


10.47216149547027 


3007r 


10.47216132022817 


10.47216132022820 


10.47216132022773 


BSOtt 


10.47216138600882 


10.47216138600883 


10.47216138600868 



Table 2: Comparison of the value of Im(0(25)) for different values of parameters 
N and d. 

evaluated at the optimum a* for higher computer precisions. The digits in 
bold correspond to correct digits of the Stokes constant. We also note that the 
numerics suggest that Qq is pure imaginary which agrees with our prediction. 

Finally, let us mention that in the process of computing the Stokes constant 
we have made several choices for the parameters. Namely, the number of terms 
N used to compute F^v and the parameter d which were used in computing the 
initial conditions of step (ii) of the numerical scheme. In fact the results are 
independent of these particular choices and Table [5] demonstrates the robustness 
of the numerical method. 

3 High precision computations of an asymptotic 
expansion for the homochnic invariant 

In this section we present a numerical method for the computation of the homo- 
clinic invariant as defined in ([25]) for the Swift-Hohenberg equation with k — 2 
and e < 0. Moreover we investigate from a numerical point of view the valid- 
ity of the asymptotic expansion (1321) for the homoclinic invariant. This section 
follows the ideas of [8] originally developed for the study of exponentially small 
phenomena for area-preserving maps. 

In order to compute the homoclinic invariant ([9]) we need to compute two 
tangent vectors at the symmetric homoclinic point r|(0,0). Using the fact 
that the system is reversible we can obtain the stable tangent vector d^pT^ 
by applying the reversor to the unstable tangent vector dtpV^. The unstable 
tangent vector d^V^ lives in the tangent plane of the unstable manifold at the 
symmetric homoclinic orbit. Thus an easy way to compute this tangent vector 
is to approximate the primary homoclinic orbit near the equilibrium point by 
the following expansion, 

N / k \ 

^lNi'P,z) = X^e'^M Cfc(e) +^afej(e)cosOV) -l-bfej(e)sinOV) (37) 
fc=i \ j>i I 
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and then use the variational equations, 

v' = DXh, (x)v ^ ' 

to transport the tangent vector d^T^ along the primary homoclinic orbit until 
it hits the symmetric plane Fix(S') defined by {(72 = 0,pi = 0}. Let us present 
the details of the method. 



3.1 A method for the computation of the homochnic in- 
variant 

(i) The first step is to determine the coefficients of ([57]) . To that end we take 
a new expansion, 

N / k 
un{(P, 2;) = ^ e''^ Cfc(e) + ^ afcj(e) cos^V) + 6fcj(e) sm{jip) 
k=i \ j>i 

and substitute into the equation, 

{{a,d^ + Mzf + lfu^eu + 2u^ - (39) 

and collect the terms of the same order in e^'^ . In this way it is possible 
to determine coefficients Cfe, a^j and b^j. It is not difficult to see that the 
coefficients ai^i and bi^i satisfy no relations and that all other coefficients 
depend from these two. So we define, 

ai,i = rocos('0o) and 61,1 = rosin('0o) 
Now recall that the first component of solves equation ([59)1 and due to 



the asymptotic behavior (p2)) we conclude that for z << and 5 « \ li 
is approximately, 

^ cos((^) + ^ f 1 + ^) sin((^)) + 0(e2^) (40) 

where e = —AS^. Next we "match" the leading order of U]^{(j>, s) with the 



(41) 



expression (|40p and conclude that ipo and tq must satisfy, 
ipo = arctan 

ro = 




Taking into account ^ we reconstruct F"^ from ujy and due to the 
"matching" (|4ip we have, 

T'^{t,t) ^Tlj^{t,t), as S^O. 

That is, for small values of (5, the expansion F^jy provides a good approx- 
imation of the primary homoclinic orbit near the equilibrium point. 
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(ii) The second step is to improve the accuracy of the approximation of the 
symmetric homochnic point, provided by F"^. Given small 5 and suffi- 
ciently large Tq > we want to determine (T, ^) such that, 

x' = XhM), x(0; V) = Tlj^{-a,To,-PeTo;i^) 

subject to, 

x(T; V) e Fix(5) (42) 

This problem can be solved using Newton method. Starting from (Tq, ipo) 
we obtain a sequence of points {Ti,ipi), 

WiJ W \^iT^;i^^) yPim;A)J ^ ' 

that converges to a limit (T*,?/',) such that x.{T^,;'ip*) G Fix(S'), provided 
(TqjV'o) is sufficiently close to (T*,-)/)*) (see pi). The derivatives in P5]) 
can be computed using the variational equations along the orbit x(i;^). 
Later we will see that the formulae (|41l) provide sufficiently accurate initial 
guesses yielding the convergence of the Newton method. 

(iii) Having obtained in the previous step an accurate approximation of the 
symmetric homoclinic point, the last step is to integrate numerically the 
system, 

x' - XhA^), x(0;i^) = r^_^(-aeTo,-/3,To; V*) 

v' = DXhA^}^, v(0;7/;) = a,d^Tl^{-a,To,-p,To;i;,) 

and evaluate the homoclinic invariant, 

Lo = ri(v(r*, -0*), s'(v(T*, ■0*))) 

3.2 Numerical results 

We have considered a finite set I consisting of points in the interval e € 
[~jqi~j^qq] and computed the homoclinic invariant for those points using the 
method previously described. For all points in I the magnitude of homoclinic 
invariant ranges from 10~^ to 10~*^. Thus, in all numerical integrations we 
have used a high order Taylor method which allows to perform the numerical 
integration with very high precision. We have computed the coefficients of the 
expansion ((37|) up to = 5 and for each e G I we have chosen Tq sufficiently 
large so that j^{—aeTo, —/3^Tq) approximates the unstable manifold within 
the required precision. The initial point (Tq, tpo) used in Newton method proved 
to be very close to (T*,V'*) and its relative error can be observed in Figure [5l 
After computing the homoclinic invariant we have normalized it using the for- 
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-0,1 -0,08 -0,06 -0,04 -0,02 



Figure 5: Relative error of (To, ijjo) depending on e G X 










1^2 


1S3 


LJ4 


CJ5 


5 


10.47216195694 


8.979943127 


- 42.60110 


152.88 


- 774.4 


3.8x10-' 


6 


10.472161956944 


8.979943127 


- 42.601100 


152.888 


- 774.2 


3.8x10^ 


7 


10.4721619569443 


8.9799431275 


- 42.6011004 


152.887 


- 774.40 


3.80x10^ 


8 


10.47216195694439 


8.97994312752 


- 42.60110043 


152.88795 


- 774.39 


3.814x10^ 


9 


10.472161956944398 


8.9799431275209 


- 42.601100432 


152.88795 


- 774.394 


3.813x10^ 


10 


10.4721619569443983 


8.9799431275210 


- 42.601100432 


152.887958 


- 774.3944 


3.8138x10^ 


11 


10.4721619569443983 


8.9799431275210 


- 42.601100432 


152.887958 


- 774.3944 


3.813x10^ 


12 


10.4721619569443983 


8.9799431275210 


- 42.6011004327 


152.887958 


- 774.3944 


3.813x10^ 



Table 3: Coefficients of the estimated polynomials for differents subsets of V 
and different degrees 



mula, 



The behaviour of the function w(e) can be observed in Figure [6] It possible 
to see that it is approaching the value of the Stokes constant computed in 
the previous section. Moreover, it is aproaching this value in a linear fashion, 
supporting the validity of the asymptotic formula Taking into account 

the asymptotic expansion for w(e) we investigate the validity of the following 
asymptotic expansion for cD(e), 

cS(e)x^Wfce'= (44) 

fc>0 

To that end, we have taken 14 points evenly spaced in the interval [—2.7 x 
lO^-^,— 1.4 X 10"'^] and computed the corresponding normalized homoclinic in- 
variant with more than 40 correct digits. Let us denote this set of homoclinic 
invariants by V . Then, in order to get the first few coefficients of the asymptotic 
expansion (j44l) we have fitted a partial sum of the asymptotic expansion to the 
points of V ■ Here we have used as many points as the number of unknown 
coefficients. Moreover, following [S] we have performed the following tests to 
evaluate the validity of the asymptotic expansion: 

(i) Interpolating different partial sums to different subsets of V should give 
essentially the same results for the coefficients. 

(ii) The constant term of the interpolating polynomial should coincide with 
the value of the Stokes constant computed in the previous section. 

(iii) The interpolating polynomial should reasonably approximate d)(e) outside 
the interval [—2.7 x 10"'^,— 1.4 x 10~^], in the sense that it agrees with 
the main property of an aymptotic expansion: 



n-l 

w(e) - ^aifce*^ 

fc>0 



<Ce", Vee [eo,0) 



for some C > and eo < 0. 
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0,01- 



-0,1 -0,08 -0,06 -0,04 -0,02 

Figure 7: Relative error of the asymptotic expansion of w(e). 

For the first test we have considered ah possible subsets of V having only 6 
consecutive elements and interpolated these data by polynomials of degree 5. 
Then for each coefficient, we extracted the part of the number which is equal 
to all polynomials. We have repeated this process for polynomials of degree 6 
up to degree 12. The results are summarized in Table [3l where it is possible 
to see that there is a good agreement between the coefficients of the different 
interpolating polynomials of different subsets of P. We can also infer from Table 
[3] that the results are numerically stable. Thus, we have the following estimates 
for the first 6 coefficients of (j44|) : 

ujn = 10.4721619569443983 ... wi = 8.9799431275210 ... W2 = -42.601100432 . 
W3 = 152.887958 .. . W4 = -774.3944 .. . 3.813 ... x 10^ 

Furthermore, it is clear that the coefficient wq coincides (up to 18 digits) 
with the value of the Stokes constant which we recall, 

IGol = 10.47216195694439835828552143203190. . . 

Moreover, in Figure[7]we see that the relative error of the asymptotic expansion 
does not exceed 0.06 in the hole interval [— ]^,0]. Thus, our numerical results 
provide a satisfactory numerical evidence that supports the correctness of the 
asymptotic expansion (|32[) . 
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A Transformation of GSHE to the normal form 



In order to normalize He up to order 5, we have used the method of Lie series 
to determine Hamiltonians Fj, i = 0, ...,4 which generate the following near 
identity canonical map, 

where 

V 32 ^^^^ + 32 + gPiP2 

7 95 17 5 

216 216 72 36 

517 2 29 \ 3 / 217 2 17 \ , 



20736 512y^^ V 20736 512 
/ 2327 2 31 \ 2 / 19 2027 2 \ 3 
(,20736 "^ -512j«^^^^^+ (,-512 + 20736 " ' ^^^^ + 
5 7 2\ s 19 2 3 , . 

128 + 192 " r^^^+ 5^" -^)piP2^ 



/ 143 ^ 2 167 r- ,\ 
^^ = H"n52"^^^ ^^-T728 "^^^ j 

Ft = --^ V2k (37 «2 - 27) P2^ - V^'i (-45 + 52 k^) 
izio d4o 

^\/2k (-27 + 34^2)^12^2^ 



(45) 



Using an algebraic manipulator it is not difBcult to see that ^5 transforms 
into the desired form. 
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